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Abstract. Many scientific deliverables of the next generation low frequency 
radio telescopes require high dynamic range imaging. Next generation telescopes 
under construction indeed promise at least a ten-fold increase in the sensitivity 
compared with existing telescopes. The projected achievable RMS noise in the 
images from these telescopes is in the range of 1-10/iJy/beam corresponding to 
typical imaging dynamic ranges of 10^~^. High imaging dynamic range require 
removal of systematic errors to high accuracy and for long integration inter- 
vals. In general, many source of errors are directionally dependent and unless 
corrected for, will be a limiting factor for the imaging dynamic range of these 
next generation telescopes. This requires development of new algorithms and 
software for calibration and imaging which can correct for such direction and 
time dependent errors. In this paper, I discuss the resulting algorithmic and 
computing challenges and the recent progress made towards addressing these 
challenges. 



1. Introduction 

Aperture synthesis array telescopes combine signals from a number of antenna 
pairs to sample the coherence function (visibility function) in the radiation far 
field. The angular resolution is inversely proportional to the largest separation 
between the antennas (baseline) compared to the wavelength of observation. 
The sensitivity is proportional to the total collecting area, square root of the 
bandwidth of observation and total integration time and inversely proportional 
to the effective system temperature. 

Next generation telescopes radio telescopes, some under construction, promise 
10-100 times improvement in resolution and sensitivity. For a number of reasons 
ranging from engineering challenges and cost considerations to sky background 
emission, it is hard to lower the system temperature by a few order of mag- 
nitude to improve the telescope sensitivity by similar order. As a result, all 
next generation telescopes use larger number of antenna elements to increase 
the collecting area, wide-band receivers and long integrations in time to achieve 



the higher sensitivity. To mitigate bandwidth smearing (Thompson et al. 2001 1, 
effects of narrow band radio frequency interference (RFI) and for scientific rea- 
sons, the observed band is split into a number of narrower frequency channels. 
Snapshot data rate is proportional to the product of the square of the number of 
antenna elements and number of frequency channels. This fact, combined with 
long integrations in time implies that the projected sensitivity improvements of 
the next generation telescopes will come at the cost of lO'^"^ times increase in 
the data volume over existing telescopes. 
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An underlying assumption in the sensitivity calculations is that the random 
noise in the observations has no systematic component and that for a given sys- 
tem temperature, the signal to noise ratio (SNR) is proportional to \/ AT Au 
where AT is the total integration in time and Azv is the total bandwidth of 
observation. However, the observed data is inevitably corrupted by a number 
of instrumental and ionospheric/atmospheric effects. Furthermore, these effects 
are not the same across the field of view (i.e., these effects are direction de- 
pendent (DD)). This makes the "noise" in the observations non-random which 
does not necessarily reduce with integration in time and/or frequency. Further- 
more, low frequency sky is also brighter and more complex. As a result, the 
projected image plane RMS noise of l-10//Jy/beam translates to an imaging 
dynamic range requirement of 10^~^. The imaging dynamic range limit due to 
deconvolution errors for complex fields with compact and extended emission is 
significantly higher than this. 

The next generation post processing software therefore needs to correct for 
direction dependent effects more accurately, over larger parameter space (time, 
frequency and polarization) using 2-4 orders of magnitude larger data volume 
as well as image complex sky emission with high fidelity to achieve the scientific 
goals. An obvious conclusion is that we necessarily need significant research in 
the area of post processing techniques for imaging and calibration and develop 
algorithms which are more accurate, account for direction dependent effects and 
can deal with large data volumes efficiently. In this paper, I review the recent 
progress in the development of new imaging and calibration algorithms relevant 
for high dynamic range imaging at low radio frequencies (< 2 GHz). A more 
complete theoretical background can be found in the recent paper by |Rau et ah 
(20091). 



2. The Measurement Equation 



Using the theoretical formulation by Hamaker et al. 



measurements from a single baseline can be described 
ment Equation 



(19961, full polarimetric 



3y the following Measure- 



V^^'{iy,t) = Jij{u,t)Wij{u,t) I E^j{s,u,t)r{s,u,t)e'^^^-'ds 



(1) 



where 1^9'"* is the observed visibility samples measured by the pair of antennas 



designated by the subscript i and j, separated by the vector bij and weighted by 
the measurement weights Wij. Jij is the complex direction independent gain, Eij 
is the direction dependent gain as a function of the direction s, frequency and 
time t and / is the image vector. The vectors V and / are full polarization vectors 
in the data and image domain respectively. Jij and Eij can be expressed as an 



Ji{iy, t)®J*Aiy, t) 



outer product of two 2x2 antenna based Jones matrices as Jij 
and Eij = Ei{s, u, t)(^E*{s, v, t). Ji and Ei describe the full polarization response 
of the individual antennas in the feed polarization bases. An appropriate unitary 
transform can be app l ied to convert th e above equation to Stokes bases (see 
Hamaker et ah] ([l996l) ; |Rau et al.| (|2009| for details). 



For wide-band observations, the sky emission also changes as a function of 
frequency and is potentially differently for different directions. Assuming time 
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invariance, this frequency dependence can be expressed as 



(2) 



where 1° is the image at the reference frequency Uq and a is spectral index which 
varies across the field of view and the frequency band. 1°, Ji, Ei and a repre- 
sent the explicit unknowns in Equation [T] The process of calibration estimates 
Ji and Ei while the process of imaging estimates 1° and a{s,v). Note that 
while the observed visibilities can be corrected for the effects of JjS by dividing 
Equation [T] by Jjj, the same is not true for DD terms. Correction for the effects 
of Eij can only be done as part of the imaging process. This makes solving and 
correcting for DD errors more difficult and consequently conventional calibration 
accounts for only direction independent corruptions. This has been sufficient for 
the existing telescopes. This however must change to achieve imaging dynamic 
ranges consistent with the thermal noise limit of the next generation instruments 
(Bhatnagar et al.||2004l 120061 iRau et al.||2006p. 



3. Parametrization of the Measurement Equation 

In order to make an image free of the effects of the corruptions, V^^^ needs to 
be corrected for the effects of JjS and EiS. Conventional techniques typically 
parametrize Jij as three separate terms to represent time, frequency and polar- 
ization dependencies. These are assumed to be orthogonal and therefore solved 
independently in the process of time, bandpass and polarization calibration. For 
narrow band observations (less than 10% fractional bandwidth), a is assumed 
to be small or zero and 1° is parametrized as a value per pixel, each pixel being 
treated as independent degree of freedom (DoF). Ei is ignored in calibration 
and during image deconvolution and if required, corrections for it are made 
post-deconvolution (e.g. post-deconvolution correction for the antenna response 
as a function of direction) . 

For the next generation telescopes however, more sophisticated parametriza- 
tion is required. The DD terms need to be parametrized to model the instru- 
mental and ionospheric DD effects. Variations of a across observing frequency 
band and across the field of view (FoV) needs to be parametrized to model the 
spectral index of the sources. Sky is brighter and more complex at low frequen- 
cies and most fields have sources with extended emission. 1° therefore also needs 
to be parametrized to better represent extended emission. 



4. The W-term 



Uijl+Vijm+Wij 



The exponent in Equation jlj can be expanded as bij - 

with the usual meaning for the symbols ( Thompson et aL]|2001 1. When the field 
of view is large (as is typically true at low frequencies), the integral cannot be 
reduced to a Fourier transform and use of the 2D FFT algorithm for computing 
efficiency leads to significant distortions away from the phase center. 

Faceting algorithms using faceting in the uv-domain (as against image- 
plane faceting) (see Sault et al. (19991 for an expression for faceting in the uv- 
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Figure 1. Figure shows the performance of imaging algorithms to c;orrect 
for the effects of the w-term. Image on the left was made using the uv-faceting 
algorithm. Image on the right was made using the w-projection algorithm. 
Compact sources well away from the center of these images arc undistortcd. 
The RMS noise in the two images is the same. Residual errors are more 
systematic for the facet based algorithms. 
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Figure 2. Image showing wide- field P-band imaging in the Galactic plane 
using the W-Projection algorithm. Extended emission as well as compact 
sources away from the phase center show no distortions due to the w-term. 



domain) produce undistorted single-plane images of the sky with an approximate 

space-invariant PSF. This has several practical advantages during dcconvolution, 
particularly for the deconvolution of extended emission. Algorithms using this 
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approach exist in the CASA pac kagj^ Computing load is the s ame as that for 
image-plane faceting algorithms ( Cornwell Perley||1992 1999). 



4.1. The W-Projection Algorithm 

Absorbing the third term in the expression for bij ■ s into Eij in Equation [l] the 

observed visibilities can be expressed as V^''' = [FT (Eij)] where V are the 
visibilities corresponding to the tangent plane, FT represents the Fourier trans- 
form and denotes the convolution operations respectively. The W-Projection 
algorithm ( [Cornwell et al. 20081 exploits this to correct for the w-term in the 



gridding/de-gridding operation during imaging. While theoretically this algo- 
rithm can be shown to be faster by up to 50 times, in practice it has been shown 
to be up to an order of magnitude faster compared to faceting algorithms. W- 
Projection algorithm also produces a single-plane image, making it easier to 
combine with other techniques for dealing with extended emission across the 
field of view as well as correcting for other DD effects. 

5. Ionospheric corruptions 

Corruptions due to ionosphere is one of the limiting problems in high sensitivity 
high resolution imaging at low frequencies. Its effect is that the phase across the 
antenna aperture is not constant and potentially different for each antenna in 
the array (it is a direction dependent effect - often referred to as "non-isoplanatic 
ionosphere" in the literature). 

5.1. Field base calibration 



The field-based calibration technique ( Cotton et al.|[2004 | measures the shift of 



compact sources throughout the FoV to estimate local ionospheric phase gradi- 
ents. A polynomial fit to the estimated phases is then used to apply corrections 
to the rest of the field. For small baselines (< 2 — 3Km) this has been shown to 
improve the imaging performances for some fields. This technique however can 
be computationally prohibitive for large FoV with complex emission and does 
not extend to cases where the ionospheric refractive effects are significant. 

5.2. Peeling 

This technique estimates a complex gain towards a number of sources across the 
FoV for which good models are known apriori (either from earlier imaging and 
calibration runs or from external sources). The solved gains are then used to 
remove the contribution of the sky emission in the vicinity from the data as: 

yCorrected ^ ^Obs _ ^ jPy.f " (3) 

k 

where the superscript k denotes all sources in the region where the peeling so- 
lutions Jj^ apply- In the iterative form of Peeling, the corrected visibilities are 



^CASA Home Page at 



http : //casa . nrao . edu 
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then used to apply this technique iteratively to the strongest sources in the 
residual image. This has been shown to work well for simple fields (dominated 
by compact strong sources) and for relatively small data volume. Various vari- 
ants of this techniq ue are being currently tested (Nijboer Sz Noordam 2007[ 
Mitchell et al.|[2"008 ) to determine its numerical and computational performance 



for complex fields with extended emission and with large data volumes. 



6. Effects of antenna Primary Beam 

For narrow band observations, time dependent DD gains are dominantly due 
to time varying antenna primary beams (PB). While the antenna forward gain 
is clearly direction dependent, its time variation is due to a number of reasons 
(rotation of the rotationally asymmetric PBs with Parallactic Angle for Az-El 
mount antennas, antenna pointing errors, geometrical distortions of the antenna 
with elevation, etc.). For wideband observations, the shape of the PB varies 
across the band and sources well within the PB main-lobe at the lower frequency 
end of the band may appear in the first sidelobe at the higher frequency end (e.g. 
with a bandwidth ratio of 2:1). For sources in the first sidelobe of the antenna 
power pattern, the time varying gain due to the rotation of the PBs will be 



even stronger (Bhatnagar et al. 2006, 20081. Furthermore, the time, frequency 



and direction dependence of aperture array station power pattern is expected 
to be worse compared to filled aperture antennas. Algorithms to correct for PB 
effects are therefore crucial for the scientific deliverables of the next generation 
instruments. 

Image-plane based PB correction by direct evaluation of the integral in 
Equation [T] is possible. To reduce the resulting prohibitive run-time comput- 
ing cost for realistic image complexity and data volumes, a FFT based reverse 



( 2008 ) . However this requires mak- 



ler the sky emission or the antenna 



transform has been used by Uson & Cotton 
ing assumptions about the variability of eit 
power pattern. 

6.1. The A-Projection Algorithm 

Eij represents the effects of antenna primary beams in Equation [Tj The A- 
Pro jection algorithm ( Bhatnagar et al.]|2008 1 uses a model for antenna aperture 



illumination and the approximate unitary nature of the resulting operator to 
correct for effects of PB as part of image deconvolution iterations. This algo- 
rithm can naturally deal with non-identical antenna PBs, is straight forward to 
integrate with other advanced imaging and calibration algorithms and is com- 
putationally efficient. An estimate of the antenna aperture illumination pattern 
is however required - which can be measured (antenna holography) or modeled. 
See Fig. |3]for an example of application of this algorithm for imaging at L-Band 
using the VLA. 



7. The Sky Model 

The sky model / in Equation [l] is computed using image deconvolution algo- 
rithms. Most conventional algorithms treat each pixel with significant emission 
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Figure 3. Figure shows the results of the application of the A-Projection 
algorithm for VLA L-Band imaging. Top panel shows Stokes-I images made 
using conventional (left) and the A-Projection algorithm (right). Bottom 
panel shows Stokes-V images. Stokes-V imaging with the VLA suffers from 
strong and time varying instrumental effects which are completely corrected 
in the image in bottom right panel. 



in the image as an independent DoF (e.g., the Clean ( Hogbom||1974 | and MEM 
( Cornwell Evans|1985 1 algorithms and their vari ants). Such a parametriz ation 
of the sky is non-optimal for extended emission (Bhatnagar &; Cornwell 2004 1 
and suffers from the problem of pixel quantization errors ( [Voronkov Hi Wieringa 



2004 



prob 



Cotton Sz Uson|2008 1 . For complex fields with strong emission, both these 



ems limit the imaging dynamic range well above the instrumental limit of 



the next generation telescopes. 



7.1. Scale-sensitive modeling: The MS- and Asp-Clean algorithms 

The MS-Clean algorithm ( CornweTIl|2008 1 models the sky as a collection of com- 
ponents with pre-computed set of scale sizes. Depending upon the user defined 
choice of scales, extend emission as well as compact emission can be better mod- 
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Figure 4. Figure showing the performance of MS-Clean and Asp-Clean al- 
gorithms for the deconvolution of complex extended emission. Image on the 
left is the model image used for the tests. Images in the middle and on 
the right show the residuals from MS-Clean and Asp-Clean algorithms re- 
spectively. Spatial correlation scale is significantly reduced in both residuals 
compared to residuals from the Clean algorithm (not shown) . Residuals from 
Asp-Clean algorithm are more noise like compared to MS-Clean algorithm. 



eled with a significantly smaller number of components (DoF). Memory require- 
ments and computing load is higher compared to conventional algorithms and 
the coupling between fixed set of user defined scales is ignored (i.e., it effectively 
ignores the fact that the parameter space is non-orthogonal). 

The Asp-Clean algorithm (Bhatnagar &; Cornwell 2004) adaptively deter- 
mines the local scale as well as position of the components in the image. This 
effectively mitigates the problem of pixel quantization and accounts for cou- 
pling between various scales (i.e., recognize the fact that parameter space is 
non-orthogonal). Although for the same number of components the comput- 
ing requirements are 2-3 times higher compared to MS-Clean, the number of 
components required is significantly smaller for same image complexity. 



7.2. Wide-band modeling: The MS-MFS algorithm 



Ignoring the frequency dependence of the sky in wide-band o bservation can 
limit the imaging dynamic range to ~ 1 : 10^ (Rau et al. 20061. Hence, apart 
from scale-sensitive modeling of the sky, modeling the frequency dependence 
of the sky (Equation l2| during imaging is required for wide-band observations. 
The MS-MFS algorithm pill|| 2009] ) uses the MS-Clean approach to model ex- 
tended emission and models the frequency dependence using a Taylor expansion 
of Equation [2] about the reference frequency and solving for the coefficients of 
the series. The A-Projection algorithm is used to correct for the PB frequency 
dependence in combination with MS-MFS to make S tokes-I, spectral index and 
spectral index variation images of the sky ( Rau| 20"09 |. The combined algorithm 
is being currently tested using wide-band observations of fields with strong ex- 
tended emission. 
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8. Solvers for direction dependent effects 



For high resolution high dynamic range imaging, it is virtually impossible to 
measure the DD terms in Equation [T] to the required accuracy prior to imaging. 
Algorithms to model and solve for the DD effects are therefore required. 

Peeling based algorithms attempt to solve for DD effects by allocating few 
DoF per direction of interest (Dol). Solutions for each Dol are either used locally 
to subtract from the data (section [5.2. ), or interpolated for other sources as well. 




Figure 5. Figure showing simulations with typical pointing errors for VLA 
antennas as a function time (red curves drawn with lines and symbols). The 
solutions for antenna pointing errors derived using the Pointing SelfCal algo- 
rithm are the over-plotted curves (blue). 



8.1. Pointing SelfCal 

Another approach, fundamentally different from Peeling, is to develop physical 
models for the various DD effects and solve for the parametric model using 
Equation [!} Projection methods to correct for known DD effects (sections 4.1 



and 6.1.[ ) can be easily used to implement solvers which are also computationally 



efficient. This approach fundamentally mitigates the problem of the proliferation 
of DoFs inherent in the peeling approach. The Pointing SelfCal algorithms 
dBhatnagar et al.|2004p is an example of use of this approach to solve for antenna 



pointing errors. To correct for the pointing errors, the solved pointing errors are 
included as part of the model for the antenna aperture illumination and used 
in A-Projection algorithm during imaging. Fig. |5] shows results from tests using 
simulated data. Further work on this using real data is currently in progress. 
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